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ABSTRACT 

We study the escape rate of stars, N, from clusters with different radii on circular 
orbits in a tidal field using analytical predictions and direct iV-body simulations. We 
find that N depends on the ratio £H = fh/Vj, where rh is the half- mass radius and 
rj the radius of the zero-velocity surface around the cluster. For > 0.05, the "tidal 
regime" , there is almost no dependence of N on 9\. To first order this is because the 
fraction of escapers per half-mass relaxation time, t r h, scales approximately as 
which cancels out the rf/ 2 term in t r h. For 9i < 0.05, the "isolated regime", N scales 

as 9l~ 3 / 2 . The dissolution time-scale, idis, falls in three regimes. Clusters that start 
with their initial 91, 91^, in the tidal regime dissolve completely in this regime and 
their tdis is, therefore, insensitive to the initial r^. Our model predicts that 9^ has 
to be 10~ 20 — 10~ 10 for clusters to dissolve completely in the isolated regime. This 
means that realistic clusters that start with 9^ < 0.05 always expand to the tidal 
regime before final dissolution. Their idis has a shallower dependence on 91^ than what 
would be expected when idis is a constant times i r h- For realistic values of 9^, the 
lifetime varies by less than a factor of 1.5 due to changes in !EHj. This implies that the 
"survival" or "vital" diagram for globular clusters should allow for more small clusters 
to survive. We note that with our result it is impossible to explain the universal peaked 
mass function of globular cluster systems by dynamical evolution from a power-law 
initial mass function, since the peak will be at lower masses in the outer parts of 
galaxies. Our results finally show that in the tidal regime <dis scales as N - 65 /lu, with 
u> the angular frequency of the cluster in the host galaxy. 

Key words: stellar dynamics - methods: TV-body simulations - globular clusters: 
general - galaxies: star clusters 



1 INTRODUCTION 



Stars escape from clusters due to internal two and three 
body encounters in which stars get accelerated to veloc- 
ities higher than the escape velocity. The resulting dis- 
solution time-scale, tdis> depends on the number of stars, 
N, and the escape rate, N, as t^is = —N/N. For a con- 
stant N the instantaneous value of idis is the remain- 
ing time to total dissolution. For clusters in isolation 
tdia scales linearly with the half-mas s relaxation time, i r h 
l|Ambartsumianl Il938l ; ISpitzerl Il940l ). The presence of a 
tidal fie l d spe e ds up N by roughly an order of magnitude 
ilHenonl Il96ll : ISpitzer fc Chevalier! 1 19731 : iGiersz fc Heggiel 
1 19971 ). When treating the tidal field as a radial cut-off, 
as is often done in Fokker Planck and A-body simula- 
tions dChernoff fc Weinberg||l990l: [Gnedin fc Ostrikerlll997l ; 
Taka hashi fc Portegies Zwart|[2000l ) idis also scales with t r h 
|Baumgardtll200ll , hereafter B01). 



iFukushige fc Heggiel (|2000h demonstrated that for a re- 
alistic tidal field it is important to consider the finite time 
it takes stars to escape through one of the Lagrange points. 
B01 showed that then tdis oc t^ 4 . This scaling was also found 
for more realistic A^-body simulations that include a stel- 
lar mass function and stellar evolu tion (|Vesperini fc Heggiel 
1 19971 : iBaumgardt fc Makindl2003l . hereafter BM03). BM03 
only considered the dependence of t T h on cluster mass, M, 
since in their runs the initial half-mass radius, r^, and the 
initial tidal radius are linked, so oc M 1 ^ 3 . 

Observations of young (extra-galactic) clusters show 
that the scaling of rh with cluster mass, M, and galactocen- 
tric distance, Ra, is considerably shallower: rh oc 
i|Zepf et al l 1199a iLarsenl 120041 : IScheepmaker et all I2007D 
than the Roche- lobe filling relation (r\ oc M 1 / 3 R^ 3 ), im- 
plying that massive clusters and clusters at large Rg are 
initially under-filling their Roche-lobe. 

The relation betwee n tdis and the median cluster ra- 
dius, f, was modelled by IWielenl { 197ll ) for clusters up to 
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TV = 250. He found that t d is (in Myrs) scales with f 3/2 for 
r < 0.01 rj, with rj the Jabobi radius, being the radius of 
the zero- velocity surface around the cluster imposed by the 
tidal field. This is because for these clusters tdis oc i r h. For 
f > 0.05 r j this scaling is not followed any more and tdis is 
shorter. Large clusters are even more vulnerable to disrup- 
tion when the effect of passing molecular clo uds is included 
llKindl 19581; I Wielenll 19851; iGieles et al.ll2006l ). 

iTanikawa fc Fukushigel (|2005r ) modelled the evolution 
of clusters of larger TV that are initially Roche-lobe under- 
filling using collisionless TV-body simulations. They confirm 
that for Roche- lobe filling clusters tdis oc tf h , with x — 3/4 
as was found by B01, and x somewhat larger for Roche- 
lobe under-filling clusters. They do not discuss the effect of 
radius on tdis- More realistic TV-body simulations, including 
various stellar initial mass function and stellar evolut ion, of 
Roche -lobe under-filling clusters were considered by lEngld 
i| 19991 ). She concludes that Roche-lobe under-filling clusters 
survive longer. However, this is for a fixed fh,i and different 
rj. Since her simulations included stellar evolution it is not 
possible to scale these results to the same rj. 

Predictions for the survival prob ability of globular 
cluster, such as the su rvival triangle (jFall fc Reesl 1 19771 ; 
iGnedin fc Ostrikeij[l997l ). rely on the assumption that t d i s is 
a constant times i r h. This is also an important assumption 
in a recent attempt to explain the shape and the univer- 
sality of the turn-over of the globular cluster mass function 
i|McLaughlin fc FaUll200Sh . To be able to judge the applica- 
bility of such models it is of importance that the relation 
between tdis and r\ is better understood for clusters with 
larger TV. 

The interplay between internal relaxation effects and 
external tidal effects is the topic of this Letter. In § [2] we 
present a set of TV-body simulations of clusters with differ- 
ent initial Roche-lobe filling factors. We introduce a simple 
analytical model for TV that includes internal dynamics and 
the external tides in § [3] In § [3] we confront our model with 
the simulations and our conclusions are presented in § [5] 



2 DESCRIPTION OF THE RUNS 

We simulate the evolution of clusters containing between 
TV = 1024 and TV = 32768 particles, without primordial 
binaries and mass-loss by stellar evolution, orbiting with 
angular frequency ui in a steady point mass tidal field to 
simulate circular orbits in a galactic potential. The stellar 
masses are randomly drawn from a power-law mass function 
with index —2.35 with the maximum mass 30 times larger 
than the minimum mass. 

F or the density distribution of the clusters we use iKind 
(|l966h models, with Wo = 5. We define an initial Roche-lobe 
filling factor as 5 = r t /r j.i, with rj,i the initial Jacobi radius 
and rt the King tida l radius, tha t is, the radius where the 
stellar density of the [King] (|l966l ) model drops to zero. We 
model four values of 5, from 0.125 to 1. The Wo parameter 
and J set the ratio of the initial r^, rh,i, and rj,i, which 
we denote by 9U (= 7v/rj )t ). Clusters with TV=[1024, 2048, 
4096, 8192] are run [16, 8, 4, 2] times to reduce statistical 
variations. We also run an TV = 4096 simulation in isolation 
up to ti/2 to determine the mass loss parameters. 

All clusters are scaled to TV-body units, such that G = 



Table 1. Summary of the TV-body simulations. 



N 




^ rj,i 


^1/2 




h/2 

*rh,i 


1024 


1 


0.186 


114 


7.20 


5.8 


1024 


0.5 


0.093 


306 


6.84 


15.7 


1024 


0.25 


0.047 


768 


6.07 


39.4 


1024 


0.125 


0.023 


1718 


4.80 


88.2 


2048 


1 


0.186 


174 


11.0 


5.1 


2048 


0.5 


0.093 


478 


10.7 


13.9 


2048 


0.25 


0.047 


1184 


9.35 


34.3 


2048 


0.125 


0.023 


2768 


7.73 


80.3 


4096 


1 


0.186 


269 


17.0 


4.4 


4096 


0.5 


0.093 


765 


17.1 


12.4 


4096 


0.25 


0.047 


1833 


14.5 


29.6 


4096 


0.125 


0.023 


4235 


11.8 


68.5 


8192 


1 


0.186 


435 


27.5 


3.9 


8192 


0.5 


0.093 


1188 


26.5 


10.6 


8192 


0.25 


0.047 


2925 


23.1 


26.1 


8192 


0.125 


0.023 


6614 


18.5 


59.0 


16384 


1 


0.186 


670 


42.3 


3.3 


16384 


0.5 


0.093 


1867 


41.7 


9.1 


16384 


0.25 


0.047 


4705 


37.2 


22.9 


16384 


0.125 


0.023 


10456 


29.2 


51.0 


32768 


1 


0.186 


1062 


67.1 


2.8 


32768 


0.5 


0.093 


3049 


68.1 


8.1 


32768 


0.25 


0.047 


8011 


63.3 


21.2 


32768 


0.125 


0.023 


17281 


48.3 


45.7 



Table 2. Results of an isolated run. 



TV v t cc t rh j n rh ti/2 

4096 0.081 209 68 3.1 1.9xl0 6 



M = 1 and E = -0.25 (|Heggie fc Mathieulll98d ). Here, E is 
the total (potential and kinetic) initial energy of the cluster. 
In these units the virial radius, r v = GM/(— 4 75), equals 
unity and the crossing time at r v is 2V2. 

The half-time, Ti/2, is the when half of the initial num- 
ber of stars have become unbound, where bound is defined 
as the number of stars within rj. We multiply t\n (in TV- 
body times) by u to compare the results of different £Hi . The 
dimensionless u)t\/2 results are equivalent to physical times 
for clusters at the same Rg ■ A summary of the runs and the 
resulting t\i% and wti/ 2 values is given in TableQ] 

A ll TV-body calculations were carried out with the 
kira ^McMillan fc Hut! Il996l : iPortegies Zwart et al.l 120011 ) 
integrator on the special purpose GRAPE-6 boards 
(Mak ing et al.l [2003) of the European Southern Observatory. 



3 ANALYTICAL MODEL FOR THE ESCAPE 
RATE OF TIDALLY LIMITED CLUSTERS 

3.1 The "classical" Ansatz 

Lets first assume that a cluster, consisting of TV stars, loses 
a constant fraction £e of its stars each trh, so that we can 
write for TV (for example [Spitzer|[l987l . hereafter S87) 
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N 



(1) 



wher e t r h is conventionally expressed as (|Spitzer fc Hartl 



ATl/2 r 3/2 

0.138— — 



(2) 

VmGlnA' 

where A = 7 A and 7 = 0.11 (jGiersz fc Hegg ie 1994a), G is 
the gravitational constant and rh is the mean stellar mass. 
The crossing time at rh, t CI , is given by 



tci ~ k (gm) 



1/2 



(3) 



where A; is a constant of order unity depending on the cluster 
density profile. 

The escape energy of stars in an isolated cluster, E^° t , is 
four times the mean kinetic energy of stars in the clusters, so 
££rit = 0.8GM/rh (S87) and the escape velocity, v e , scales 
with the stellar root mean square velocity in the cluster, 
v rlnB , as v% = 4 v 2 ms . From integration over a Maxwellian ve- 
locity distribution the fraction of stars with v 2 > 4 « 2 ms can 
be determined. We refer to this escape fraction for isolated 
clusters as £ c o and S87 showed that £ c o = 0.0074. 
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Figure 1. The fraction of escapers, £ c , for different IK for clusters 
in a tidal field. The points show the result of a numerical integra- 
tion, the full line an exponential approximation and the dashed 
line a double power-law approximation. 



3.2 The escape fraction as a function of cluster 
radius 

When a cluster evolves in a tidal field the critical energy for 
escape is 



tid _ 3GM 
*>c*- 2 ■ 



(4) 



For a point-mass galaxy, rj depends on ui and M as 

>r ' M 1/3 , (5) 



rj = 



3w s 



where ui = Vg/Rg, with Vg the circular velocity. A large u) 
means a strong tidal field, which results in a small rj. 

The ratio E^ t and E% ° t gives the relative reduction of 
the escape energy due to the tidal field (following S87) 



El 



3GM /0.8GM 



2rj 



(6) 



where we have used £R = r^/rj. Note that this F is a fac- 
tor 1.5 higher than the original definition in S87, since his 
result was based on E^f t — GM/rj, which he later refines 
to equation Q. We calculate £ c as a function of SH by nu- 
merically integrating Maxwellian velocity distributions for 
different $R to determine the fraction of stars with velocities 

« 2 ^4[i-r($n)K 2 ms . 

In Fig.[T]we show that £ c increases exponentially for in- 
creasing and can be well approximated by £ c o exp(lOJK). 
For $H > 0.05, the expression for £ e scales approximately as 
?H 3//2 , which has important consequences for N (see equa- 
tions ((TJ & (J2j|). We approximate £ c by 



£e0, 



3/2 



(7) 



where fHi = 0.05 is the boundary between the "isolated 
regime" (9t < JRi) and the "tidal regime" («R > 9ti). 

Substituting equations (J2j &; (J7j| in equation fl} and 
using equation ((5} we find for N in the tidal regime 



N 



\/3gcomA 
0.138 fR 3 / 2 



(8) 



So N is independent of rh in the regime where £ c oc 
SH 3//2 . This is because for a smaller (larger) star cluster the 
shorter(longer) t r h is balanced by the lower(higher) f;Jj]. 
Equation (JSJ also shows that N depends only marginally 
on N through In A. 



3.3 Including the escape time 

iFukushige fc Heggiel (|2000D consider the time-scale of escape 
for stars in a cluster evolving in a tidal field. This time-scale 
is non-zero because stars with energies (slightly) larger than 
the escape energy still need a finite time to find one of the 
Lagrange points, where the escape energy is lowest, to leave 
the cluster. 

B01 derives an expression for N of stars in the potential 
escaper regime, that is, with energies higher than the escape 
energy, but still trapped in the potential, and shows that it 



scales as Nt, 



-3/4 -1/4 



instead of N t T ^ (equation [TJ . We 



include £ e and define N analogous to equation {TJ , as 



N 



N 



,3/4 1/4 ' 
''rh I'cr 



(9) 



1 Ivan King noticed a remarkable similarity between this result 
and equation (53) in his 1966 paper. There he shows that the 
escape rate of stars from a Roche-lobe filling cluster is indepen- 
dent of position within the cluster. This is probably because of 
the same physical reason, but he derived it in a different way. 
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With the expressions for £ e , irh and t CI we then find for TV 
in the tidal regime, including the escape time, 



TV 



v3£eo 



In A 
0.138 



TV 



V4, 



(10) 



From a comparison between equation |8j and equation (|10[) 
we see that TV becomes TV dependent when we include the 
escape time, but is still independent of 7V The dissolution 
time-scale in the tidal regime, which we define as tjj^ = 
-N/N, is then 



.4 



—) 

In A / 



3/4 



(11) 



with A = 0.138 3/4 5H? /2 fc 1/4 /(%/3Cco). For a Wo = 5 cluster 
fc = 3.85, which together with the values for SHi and £ e o 
from § EH results in A = 0.277. The term (T V/lnA) 3/4 can 



be wel l approximated by B TV , with 77 ~ 0.6 (jLamers et al.l 
2005b). The values of B and r\ depend slightly on the value 
of 7 in A = 7TV. We find the best agreement with the 
Roche-lobe filling simulations for 7 = 0.2, which results in 
B = 0.5 and 77 = 0.65, so we can write = 0.138 N 0S5 /lu. 
This sca ling of tdis with N was als o derived from obser- 
vations (iBoutloukos fc Lamersl 120031 : lLamers et al . 2005b; 
iGieles et al.ll2005l ). To compare the model to the results of 
the simulations we derive t 1 ^ 2 in the tidal regime, t^ 2 , from 
tffi!. lLamers et all (|2005al ) show that when t dis = B TV'', with 
B a constant, then t 1/2 = (B TV /rj) (1 - [1/2]''), so 



n, 



/2 



0.077 



TV 



(12) 



We assume that the effect of the escape time is the same 
for all [ft > [fti, so that we can apply equation (|12p in this 
regime. With these relations we also assume that N in the 
pre-collapse and post-collapse phase is the same. 

3.4 Dissolution in the isolated regime 

We assume that clusters with [ft < [Hi evolve in the same 
way as clusters in isolation, up to the moment that [H be- 
comes equal to [Hi. This assumption is justified by equa- 
tion ([SJ) from which we see that for a cluster with SH = SHi 
the relative contribution of the tidal field to the escape en- 
ergy is less then 10%. In the absence of primordial bina- 
ries and mass-loss by stellar evolution rh and TV remain 
roughly constant up to the moment of core collapse, t cc 
|Baumgardt et al]|2002l . hereafter B02). After t cc , isolated 
clusters evolve in a self-similar way (t r h oc t), which results 
in two fundamen tal relations for the evolution of TV and rh 
i|Goodmanlll98i S87; B02): 



N(t) = TV C , 



rh(t) = r h , C c ( — 

«:cc 



(13) 



(14) 



where rh,cc and TV CC are rh and TV at t cc . More complicated 
rela tions, assuming a non-ze ro origin for these relations, ex- 
ist (|Giersz fe H cggic 1994bT) ■ Core collapse happens after a 
multiple number of the initial t r h, t r h,i: ice = n r hi r h,i- 

From equation (|13p we find that t\i 2 for clusters evolv- 
ing completely in the isolated regime, t^° 2 , is 



100 




0.01 



0.10 



«t=(Vn)i 



Figure 2. Prediction for the dimensionless half-time, ujti/ 2 , for 
clusters in the tidal regime (equation 1 1 2 6 and part of the inter- 
mediate regime I3.5H for different N. Results of the simulations 
are shown as dots. 



/2 



0.138 x 2 1/l, n lh \ TV <R 



.3/2 



In A 



(15) 



Because clusters expand after t cc , not all clusters that start 
with [Hi < [fti will reach t^° 2 before they reach [Hi .The max- 
imum [Hi for which equation (|15|l applies is found from equa- 
tions (|14[) & (|15|) and depends on [Hi and v as 



SK2 = % ( - 



(16) 



which for [Hi = 0.05 (Fig.© and v = 0.05-0.1 (B02) results 
in [H2 — 6x 10 -5 — 6 x 10~ . If we define complete dissolution 
as the moment where only 5% of the original number of stars 
is still bound, then the corresponding value of [H2 reduces to 
~ 10~ 20 — 10~ 10 . This implies that realistic clusters never 
dissolve completely in the isolated regime. 



3.5 Combining the isolated and the tidal regime 

Clusters that start with SH2 < [fti < [Hi evolve par- 
tially in the isolated regime and partially in the tidal 
regime. Though our model allows to numerically com- 
pute ii/2 for clusters in this regime, we simply connect 
log [wt i 1 s / ° 2 (SH2)] and log [(^^(JHi)] with a straight line. The 
slope of this line, representing the [Hi dependence of wti/2, is 
log [t t i i / d 2 (^i)Ai s / 2(^2)] / log [[Hi /[Ha]. This slope is slightly 
TV-dependent, and by using equations (|12[) fc (|15fl and the 
parameters from the isolated run (Tabled we find that it de- 
creases from 0.55 to 0.35 between TV = 1024 and TV = 32768. 
For TV = 10 6 the slope would be 0.15. 

In Fig.[2]we show the u)ti/ 2 following from our model for 
different TV and a range of two orders of magnitude around 
SHi . The results of the simulations are shown as dots and are 
discussed in the next section. 
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4 COMPARISON TO TV-BODY SIMULATIONS 

The results for ivti/2 of the simulations of clusters with 
different TV and % are presented in Fig. [2] The uti/2 re ~ 
suits for 5 = 0.5 clusters are nearly the same as those for 

5 = 1. For 5 = 0.25(0.125) the u>ii/2 values are approxi- 
mately 10%(25%) shorter compared to the Roche-lobe fill- 
ing results, whereas a scaling with t r h predicts a difference 
of a factor of 4 1 ' 5 (8 1 ' 5 ) ~ 8(23). From Table [1] we see that 
clusters that start Roche-lobe under-filling, have evolved for 
a much larger number of relaxation times than Roche-lobe 
filling clusters by the time they reach ti/2- 



5 CONCLUSIONS 

The dissolution time-scale of clusters evolving in a tidal field, 
in dimensionless units (uitdis) or in physical units, is almost 
independent of the initial half-mass radius, rh,i, when rh,i 
relative to the initial Jacobi (or tidal) radius, rj4, is larger 
than SKi(= ni.i/Vj.i) > 0.05. For clusters that start with 
5Hi < 0.05, t dis scales mildly with between SK°' 55 and 
«°- 35 for the range of N we consider and even flatter for 
larger N. Only clusters that start with SHi < 10 -4 can lose 
half their stars before they reach the influence of the tidal 
field. We find that in the tidal regime tdis is mainly deter- 
mined by N and the angular frequency: tdis oc N°' 65 /uj, 
that is, what was also found by BM03 for Roche-lobe filling, 
mul ti-mass clusters diss o lving in tidal fields. 

iGnedin fe Ostrikerl (|l997T ) construct the "vital" dia- 
gram of globular clusters, which is a triangle in M vs. r\ 
space, outside which clusters should have been destroyed. 
There are numerous globular clusters with small Th outside 
this triangle, which the authors denote as lucky survivors. 
We show that small clusters can in fact survive. 

Models that try to explain the evolution of the glob- 
ular cluster mass function (GCMF) from an initial power- 
law to a (universal) peaked distribution by stellar dynam- 
ical processes are in difficulties, since such models will al- 
ways produce less dissolu tion in the outer parts of galaxies. 
iMcLaughlin & Fall (|2008h have recently proposed a solution 
to this problem by assuming that AM /At oc M/t^ and so 
dM/dt oc y/ph (equations QJ42] with £ c constant), that is, 
their tdis is determined by internal relaxation effects only 
and is independent of the strength of the tidal field. How- 
ever, our results show that £ e is not constant and, therefore, 
dM/dt does not scale as y/~f>h- This makes the universality 
of the GCMF a problem (again), when trying to explain this 
by dynamical evolution alone. 
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